clear all;

lambda_SOC=0.1*0.5;

spin_Delta=0.0;
spin_nvec=[0,0,1];
spin_nvec=spin_nvec/norm(spin_nvec);

sigx=[0,1;1,0];
sigy=[0,-i;i,0];
sigz=diag([1,-1]);
sig0=eye(2);

gen_SOC_term=@(bdvec) sqrt(2)*(i)*lambda_SOC*(bdvec(1)*sigx+bdvec(2)*sigy+bdvec(3)*sigz);


a1=[2.7508   -1.5882    4.4921];
a2=[0    3.1764    4.4921];
a3=[-2.7508   -1.5882    4.4921];

Rh_pos=[0.125000 0.125000 0.625000 ;
    0.625000 0.125000 0.125000 ;
    0.125000 0.625000 0.125000 ;
    0.125000 0.125000 0.125000];

Rh_pos=Rh_pos-0.125;

Rh_posC=Rh_pos;
Rh_posC(:,1)=Rh_pos(:,1)*a1(1)+Rh_pos(:,2)*a2(1)+Rh_pos(:,3)*a3(1);
Rh_posC(:,2)=Rh_pos(:,1)*a1(2)+Rh_pos(:,2)*a2(2)+Rh_pos(:,3)*a3(2);
Rh_posC(:,3)=Rh_pos(:,1)*a1(3)+Rh_pos(:,2)*a2(3)+Rh_pos(:,3)*a3(3);



Rh_posC





%%

tetra1_posC=Rh_posC
tetra1_C=mean(tetra1_posC,1)

tetra2_posC=-Rh_posC
tetra2_C=mean(tetra2_posC,1)


newx_axis=zeros(4,3);
newy_axis=zeros(4,3);
newz_axis=zeros(4,3);
for ind=1:4
    zvec=tetra1_posC(ind,:)-tetra1_C;
    newz_axis(ind,:)=-zvec/norm(zvec);
end
newz_axis

newx_axis(4,:)=[1,0,0];
newx_axis(3,:)=[-1,0,0];
newx_axis(2,:)=[0.5,0.5*sqrt(3),0];
newx_axis(1,:)=[0.5,-0.5*sqrt(3),0];


for ind=1:4
    xvec=newx_axis(ind,:);
    zvec=newz_axis(ind,:);
    newy_axis(ind,:)=cross(zvec,xvec);
end


dot(newy_axis(2,:),newz_axis(2,:))


%%

tetra1=Rh_posC;
tetra1C=mean(tetra1,1);

tetra2=Rh_posC;
tetra2(1,:)=tetra2(1,:)-a3;
tetra2(2,:)=tetra2(2,:)-a1;
tetra2(3,:)=tetra2(3,:)-a2;

tetra1
tetra1C=mean(tetra1,1)
tetra2
tetra2C=mean(tetra2,1)


%%
% 
% figure(1);
% clf;
% 
% 
% hold on;
% linecc= [.5 .5 .5];
% plot_tetra_v4(Rh_posC,0*a1,linecc,2)
% plot_tetra_v4(Rh_posC,a1,linecc,2)
% plot_tetra_v4(Rh_posC,a2,linecc,2)
% plot_tetra_v4(Rh_posC,a3,linecc,2)
% 
% 
% %plot_tetra(-Rh_posC,0*a1,'r')
% plot_tetra_v4(-Rh_posC,a1,linecc,2)
% plot_tetra_v4(-Rh_posC,a2,linecc,2)
% plot_tetra_v4(-Rh_posC,a3,linecc,2)
% 
% plot_tetra_v4(-Rh_posC,a1+a2,linecc,2)
% plot_tetra_v4(-Rh_posC,a2+a3,linecc,2)
% plot_tetra_v4(-Rh_posC,a1+a3,linecc,2)
% 
% 
% 
% hexpts1=[];
% hexpts1(1,:)=Rh_posC(3,:)+a1;
% hexpts1(2,:)=Rh_posC(2,:)+a2;
% hexpts1(3,:)=Rh_posC(1,:)+a2;
% 
% hexpts1(4,:)=Rh_posC(3,:)+a3;
% hexpts1(5,:)=Rh_posC(2,:)+a3;
% hexpts1(6,:)=Rh_posC(1,:)+a1;
% 
% 
% hexpts2=[];
% hexpts2(1,:)=Rh_posC(4,:)+a1;
% hexpts2(2,:)=Rh_posC(1,:)+a1;
% hexpts2(3,:)=Rh_posC(2,:)+a3;
% hexpts2(4,:)=Rh_posC(4,:)+a3;
% hexpts2(5,:)=Rh_posC(1,:);
% hexpts2(6,:)=Rh_posC(2,:);
% 
% 
% 
% hexpts3=[];
% hexpts3(1,:)=Rh_posC(4,:)+a2;
% hexpts3(2,:)=Rh_posC(2,:)+a2;
% hexpts3(3,:)=Rh_posC(3,:)+a1;
% hexpts3(4,:)=Rh_posC(4,:)+a1;
% hexpts3(5,:)=Rh_posC(2,:);
% hexpts3(6,:)=Rh_posC(3,:);
% 
% 
% 
% hexpts4=[];
% hexpts4(1,:)=Rh_posC(4,:)+a3;
% hexpts4(2,:)=Rh_posC(3,:)+a3;
% hexpts4(3,:)=Rh_posC(1,:)+a2;
% hexpts4(4,:)=Rh_posC(4,:)+a2;
% hexpts4(5,:)=Rh_posC(3,:);
% hexpts4(6,:)=Rh_posC(1,:);
% 
% 
% 
% sig_x=[0,1;1,0];
% sig_y=[0,-i;i,0];
% sig_z=[1,0;0,-1];
% sig_0=eye(2);
% 
% gen_spinh=@(nnvec) (sig_x*nnvec(1)+sig_y*nnvec(2)+sig_z*nnvec(3))/norm(nnvec);
% 
% eva_spinS =@(state_wf) [state_wf'*sig_x*state_wf,state_wf'*sig_y*state_wf,state_wf'*sig_z*state_wf];
% 
% 
% 
% newnvec_tetra=Rh_posC(3,:)-tetra1C;
% newnvec_tetra=newnvec_tetra/norm(newnvec_tetra);
% 
% 
% 
% 
% 
% tripts=Rh_posC(1:3,1:3);
% tripts_xaxis=newx_axis(1:3,1:3);
% 
% frame_x=newx_axis(1:3,1:3);
% frame_y=newy_axis(1:3,1:3);
% 
% 
% %quiver3(tripts(:,1),tripts(:,2),tripts(:,3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% for ind=1:3
%     plt_newD2r_orbital(tripts(ind,:),tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
% %     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),frame_x(ind,1),frame_x(ind,2),frame_x(ind,3));
% %     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),frame_y(ind,1),frame_y(ind,2),frame_y(ind,3));
% % 
% %     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),tripts_xaxis(ind,1),tripts_xaxis(ind,2),tripts_xaxis(ind,3),'LineWidth',7);
% end
% 
% for ind=1:3
%     [dot(tripts_xaxis(ind,:),frame_x(ind,:)),dot(tripts_xaxis(ind,:),frame_y(ind,:))]
% end
% 
% 
% 
% 
% 
% shiftvec=a2;
% tripts=Rh_posC([1,2,4],1:3);
% %tripts_xaxis=newx_axis([1,2,4],1:3);
% 
% newnvec=Rh_posC(3,:)-tetra1_C;
% newnvec=newnvec/norm(newnvec);
% newnvecx=[1,0,0];
% newnvecy=cross(newnvec,newnvecx);
% tripts_xaxis=zeros(3,3);
% tripts_xaxis(1,:)=0.5*newnvecx-0.5*sqrt(3)*newnvecy;
% tripts_xaxis(2,:)=0.5*newnvecx+0.5*sqrt(3)*newnvecy;
% tripts_xaxis(3,:)=[-1,0,0];
% 
% tripts_xaxis=-tripts_xaxis;
% 
% frame_x=newx_axis([1,2,4],1:3);
% frame_y=newy_axis([1,2,4],1:3);
% 
% for ind=1:3
%     plt_newD2r_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
% %     quiver3(tripts(ind,1)+shiftvec(1),tripts(ind,2)+shiftvec(2),tripts(ind,3)+shiftvec(3),frame_x(ind,1),frame_x(ind,2),frame_x(ind,3));
% %     quiver3(tripts(ind,1)+shiftvec(1),tripts(ind,2)+shiftvec(2),tripts(ind,3)+shiftvec(3),frame_y(ind,1),frame_y(ind,2),frame_y(ind,3));
% %     quiver3(tripts(ind,1)+shiftvec(1),tripts(ind,2)+shiftvec(2),tripts(ind,3)+shiftvec(3),tripts_xaxis(ind,1),tripts_xaxis(ind,2),tripts_xaxis(ind,3),'LineWidth',7);
% end
% 
% for ind=1:3
%     %[dot(tripts_xaxis(ind,:),frame_x(ind,:)),dot(tripts_xaxis(ind,:),frame_y(ind,:))]
% end
% 
% 
% %quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% 
% 
% R3mat=zeros(3);
% R3mat(3,3)=1;
% R3mat(1,1)=-0.5;
% R3mat(2,2)=-0.5;
% R3mat(1,2)=0.5*sqrt(3);
% R3mat(2,1)=-0.5*sqrt(3);
% 
% 
% tripts_xaxis=(R3mat*tripts_xaxis')';
% shiftvec=a1;
% tripts=Rh_posC([3,1,4],1:3);
% %quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% 
% frame_x=newx_axis([3,1,4],1:3);
% frame_y=newy_axis([3,1,4],1:3);
% 
% for ind=1:3
%     plt_newD2r_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
% end
% 
% 
% tripts_xaxis=(R3mat*tripts_xaxis')';
% shiftvec=a3;
% tripts=Rh_posC([2,3,4],1:3);
% 
% frame_x=newx_axis([2,3,4],1:3);
% frame_y=newy_axis([2,3,4],1:3);
% %quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% for ind=1:3
%     plt_newD2r_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
% end
% 
% 
% 
% 
% linecolor=[0.0,0.4,0.0];
% linewidth=7;
% 
% tri_pts=zeros(3,3);
% tri_pts(1,:)=Rh_posC(1,:);
% tri_pts(2,:)=Rh_posC(2,:);
% tri_pts(3,:)=Rh_posC(3,:);
% % pp=fill3(tri_pts(:,1),tri_pts(:,2),tri_pts(:,3),'b');pp.FaceAlpha=0.5;
% 
% plot_polygon_line(tri_pts,linecolor,linewidth)
% 
% tri_pts=zeros(3,3);
% tri_pts(1,:)=Rh_posC(1,:)+a1;
% tri_pts(2,:)=Rh_posC(3,:)+a1;
% tri_pts(3,:)=Rh_posC(4,:)+a1;
% % pp=fill3(tri_pts(:,1),tri_pts(:,2),tri_pts(:,3),'b');pp.FaceAlpha=0.5;
% plot_polygon_line(tri_pts,linecolor,linewidth)
% 
% 
% tri_pts=zeros(3,3);
% tri_pts(1,:)=Rh_posC(1,:)+a2;
% tri_pts(2,:)=Rh_posC(2,:)+a2;
% tri_pts(3,:)=Rh_posC(4,:)+a2;
% % pp=fill3(tri_pts(:,1),tri_pts(:,2),tri_pts(:,3),'b');pp.FaceAlpha=0.5;
% plot_polygon_line(tri_pts,linecolor,linewidth)
% 
% 
% tri_pts=zeros(3,3);
% tri_pts(1,:)=Rh_posC(2,:)+a3;
% tri_pts(2,:)=Rh_posC(3,:)+a3;
% tri_pts(3,:)=Rh_posC(4,:)+a3;
% plot_polygon_line(tri_pts,linecolor,linewidth)
% 
% 
% linecolor=[0.0,0.8,0.0];
% 
% pt1=Rh_posC(2,:);
% pt2=Rh_posC(4,:)+a1;
% line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)
% pt1=Rh_posC(3,:);
% pt2=Rh_posC(4,:)+a2;
% line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)
% 
% pt1=Rh_posC(1,:);
% pt2=Rh_posC(4,:)+a3;
% line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)
% 
% pt1=Rh_posC(3,:)+a1;
% pt2=Rh_posC(2,:)+a2;
% line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)
% 
% pt1=Rh_posC(1,:)+a2;
% pt2=Rh_posC(3,:)+a3;
% line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)
% 
% 
% pt1=Rh_posC(2,:)+a3;
% pt2=Rh_posC(1,:)+a1;
% line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)
% 
% 
% 
% hold off;
% 
% 
% box on;
% axis equal;
% 
% 
% %view([0,0,1])
% view([0,13])
% 
% set(gca,'visible','off')
% 
% % camlight('headlight')


%%

figure(1);
clf;


hold on;
linecc= [.5 .5 .5];
plot_tetra_v4(Rh_posC,0*a1,linecc,2)
%plot_tetra_v2(Rh_posC,a1,linecc,2)
%plot_tetra_v2(Rh_posC,a2,linecc,2)
%plot_tetra_v2(Rh_posC,a3,linecc,2)


linecolor=[0.0,0.4,0.0];
linewidth=7;

tri_pts=zeros(3,3);
tri_pts(1,:)=Rh_posC(1,:);
tri_pts(2,:)=Rh_posC(2,:);
tri_pts(3,:)=Rh_posC(3,:);
% pp=fill3(tri_pts(:,1),tri_pts(:,2),tri_pts(:,3),'b');pp.FaceAlpha=0.5;

plot_polygon_line(tri_pts,linecolor,linewidth)





tripts=Rh_posC(1:3,1:3);
tripts_xaxis=newx_axis(1:3,1:3);

frame_x=newx_axis(1:3,1:3);
frame_y=newy_axis(1:3,1:3);


%quiver3(tripts(:,1),tripts(:,2),tripts(:,3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
for ind=1:3
    plt_newD2r_orbital(tripts(ind,:),tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
%     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),frame_x(ind,1),frame_x(ind,2),frame_x(ind,3));
%     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),frame_y(ind,1),frame_y(ind,2),frame_y(ind,3));
% 
%     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),tripts_xaxis(ind,1),tripts_xaxis(ind,2),tripts_xaxis(ind,3),'LineWidth',7);
end




box on;
axis equal;


view([0,0,1])

set(gca,'visible','off')


material dull
camlight('headlight')



%%
% 
figure(2);
clf;


hold on;
linecc= [.5 .5 .5];
plot_tetra_v4(Rh_posC,0*a1,linecc,2)
%plot_tetra_v2(Rh_posC,a1,linecc,2)
%plot_tetra_v2(Rh_posC,a2,linecc,2)
%plot_tetra_v2(Rh_posC,a3,linecc,2)


linecolor=[0.0,0.4,0.0];
linewidth=7;

tri_pts=zeros(3,3);
tri_pts(1,:)=Rh_posC(1,:);
tri_pts(2,:)=Rh_posC(2,:);
tri_pts(3,:)=Rh_posC(3,:);
% pp=fill3(tri_pts(:,1),tri_pts(:,2),tri_pts(:,3),'b');pp.FaceAlpha=0.5;

% plot_polygon_line(tri_pts,linecolor,linewidth)




frame_x=newx_axis(1:3,1:3);
frame_y=newy_axis(1:3,1:3);

tripts=Rh_posC(1:3,1:3);
tripts_xaxis=newx_axis(1:3,1:3);


% 
% 
% tripts=Rh_posC(1:3,1:3);
% tripts_xaxis=newx_axis(1:3,1:3);
% 
% frame_x=newx_axis(1:3,1:3);
% frame_y=newy_axis(1:3,1:3);


%quiver3(tripts(:,1),tripts(:,2),tripts(:,3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
for ind=3
    plt_newD2r_orbital(tripts(ind,:),tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
%     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),frame_x(ind,1),frame_x(ind,2),frame_x(ind,3));
%     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),frame_y(ind,1),frame_y(ind,2),frame_y(ind,3));
% 
%     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),tripts_xaxis(ind,1),tripts_xaxis(ind,2),tripts_xaxis(ind,3),'LineWidth',7);
end


% 
% %quiver3(tripts(:,1),tripts(:,2),tripts(:,3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% for ind=3
%     %plt_porbital(tripts(ind,:),tripts_xaxis(ind,:));
%     plt_D1_orbital(tripts(ind,:),tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
% end
plt_newD2r_orbital(Rh_posC(4,:),-tripts_xaxis(ind,:),newx_axis(4,:),newy_axis(4,:));
%plt_porbital(Rh_posC(4,:),-tripts_xaxis(ind,:));

linecolor=[0.0,0.8,0.0];
linecolor='cyan'

pt1=Rh_posC(3,:);
pt2=Rh_posC(4,:);
line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)


-Rh_posC

mirroc_C=mean(Rh_posC,1);
mirror_norm= Rh_posC(4,:)- Rh_posC(3,:);
mirror_norm=mirror_norm/norm(mirror_norm)
mirror_n10=[1,0,0];
mirror_n20=cross(mirror_norm,mirror_n10);

mirror_n1=mirror_n10+mirror_n20;
mirror_n1=mirror_n1/norm(mirror_n1);
mirror_n2=-mirror_n10+mirror_n20;
mirror_n2=mirror_n2/norm(mirror_n2);

%plt_porbital_convertD_v2(-Rh_posC(1,:)+a2,mirror_n1,mirror_norm,mirror_n1);
%plt_porbital_convertD_v2(-Rh_posC(1,:)+a2,mirror_norm,mirror_norm,mirror_n1);



mirror_pts=zeros(4,3);
n1l=2.0;
n2l=2.5;
mirror_pts(1,:)=-n1l*mirror_n1-n2l*1*mirror_n2+mirroc_C;
mirror_pts(2,:)=n1l*mirror_n1-n2l*1*mirror_n2+mirroc_C;
mirror_pts(3,:)=n1l*mirror_n1+n2l*0.5*mirror_n2+mirroc_C;
mirror_pts(4,:)=-n1l*mirror_n1+n2l*0.5*mirror_n2+mirroc_C;

pltmirror=mirror_pts;

hold on;
pp=fill3(pltmirror(:,1),pltmirror(:,2),pltmirror(:,3),'g');
pp.FaceAlpha=0.3;





box on;
axis equal;

%view([0,0,1])
%view([35,16])
%view([-141,-21])
view([-30,35])

set(gca,'visible','off')
material dull

camlight('headlight')




%%




figure(1)


set(gca,'visible','off');
set(gcf, 'color', 'none');
set(gca, 'color', 'none');


exportgraphics(gcf,'Fig4f-D2-3site-cancel.eps','ContentType','image','BackgroundColor','none','Resolution',400)
% 'ContentType','vector'





%%



figure(2)

view([-30,35])

view([-33,37])

set(gca,'visible','off');
set(gcf, 'color', 'none');
set(gca, 'color', 'none');


exportgraphics(gcf,'Fig4f-D2-2site-cancel.eps','ContentType','image','BackgroundColor','none','Resolution',400)
% 'ContentType','vector'



